#ifndef __INITIAL_CONDITIONS_HPP__
#define __INITIAL_CONDITIONS_HPP__
#include "Solver.hpp"
#include "stencils.h"
#include "uns2d_mesh.h"
// #include "fluidProps.h"
#include <cmath>
#include <iostream>

using namespace std;

extern vector<array<double, NQ>> consv;
extern vector<array<array<double, NQ>, DIM>> fluxp;
extern vector<array<double, NQ>> fluxe;
extern vector<array<double, NQ>> res;
extern array<double, NQ> resnorm;

/**
 * @brief 初始化流场
 * 涡波干扰， x[-65,15],y[-40,40];涡核[4,0];
 *
 */
void setInitialConditions(TYPE_MESHFSM &mesh, TYPE_STENCILS &stencils,
                          double fill_time) {
  const double x_c = 4.0;
  const double y_c = 0.0;
  const double gamma = 1.4;
  const double M_v = 1.0;
  const double R = 1.0;
  const double D = 80.0;

  // Post-shock condition.
  const double rho_post = double(1.34161490);
  const double p_post = double(1.51333333) / gamma;
  const double u_post = double(-0.89444445);
  const double v_post = double(0);

  // Pre-shock condition.
  const double rho_pre = double(1);
  const double p_pre = double(1) / gamma;
  const double u_pre = double(-6) / double(5);
  const double v_pre = double(0);

  double rho, p, u, v;

  for (int in = 0; in < mesh.nodes_.size(); ++in) {
    const double x_position = mesh.node_xyz_[in][0];
    const double y_position = mesh.node_xyz_[in][1];

    if (x_position < 1e-5) {
      consv[in][0] = rho_post;
      consv[in][1] = rho_post * u_post;
      consv[in][2] = rho_post * v_post;
      consv[in][3] = p_post / (gamma - 1.0) +
                     0.5 * rho_post * (u_post * u_post + v_post * v_post);
    } else {
      const double r = sqrt((x_position - x_c) * (x_position - x_c) +
                            (y_position - y_c) * (y_position - y_c));
      if (r > 4) {
        consv[in][0] = rho_pre;
        consv[in][1] = rho_pre * u_pre;
        consv[in][2] = rho_pre * v_pre;
        consv[in][3] = p_pre / (gamma - 1.0) +
                       0.5 * rho_pre * (u_pre * u_pre + v_pre * v_pre);
      } else {
        rho = pow(1.0 - 0.5 * (gamma - 1) * M_v * M_v *
                            exp(1.0 - r * r / (R * R)),
                  1.0 / (gamma - 1.0));
        p = 1.0 / gamma *
            pow(1.0 -
                    0.5 * (gamma - 1) * M_v * M_v * exp(1.0 - r * r / (R * R)),
                gamma / (gamma - 1.0));
        u = u_pre - M_v * exp(0.5 * (1.0 - r * r / (R * R))) *
                        (y_position - y_c) / (0.0125 * D);
        v = v_pre + M_v * exp(0.5 * (1.0 - r * r / (R * R))) *
                        (x_position - x_c) / (0.0125 * D);
        consv[in][0] = rho;
        consv[in][1] = rho * u;
        consv[in][2] = rho * v;
        consv[in][3] = p / (gamma - 1.0) + 0.5 * rho * (u * u + v * v);
      }
    }
  }
}
#endif // !__INITIAL_CONDITIONS_HPP__